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Data assimilation refers to the process of obtaining an estimate of a system's state using a model 
for the system's time evolution and a time series of measurements that are possibly noisy and 
incomplete. However, for practical reasons, the high dimensionality of large spatiotemporally chaotic 
systems prevents the use of classical data assimilation techniques. Here, via numerical computations 
on the paradigmatic example of large aspect ratio Rayleigh-Benard convection, we demonstrate the 
applicability of a recently developed data assimilation method designed to circumvent this difficulty. 
In addition, we describe extensions of the algorithm for estimating unknown system parameters. 
Our results suggest the potential usefulness of our data assimilation technique to a broad class of 
situations in which there is spatiotemporally chaotic behavior. 
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Estimation of the state of an evolving dynamical sys- 
tem from measurements that are possibly noisy and in- 
complete is a prerequisite for prediction and control in 
many situations. Furthermore, scientific investigations of 
dynamical processes often require inference of an evolv- 
ing system state from data. The term 'data assimilation' 
refers to the case where a model for the time evolution 
of the system is available and is used in conjunction with 
incoming measurements to estimate the evolving system 
state. The Kalman filter [l|, Q optimally solves the data 
assimilation problem for systems with linear dynamics. 
The classical adaptation of the Kalman filter to nonlin- 
ear systems is called the extended Kalman filter 3] . How- 
ever, it requires inversion of N x N matrices, where N 
is the number of model variables [4J. As a consequence, 
the application of such a technique to large, dynamically 
high-dimensional spatiotemporally chaotic systems is in- 
feasible because no existing computers are large or fast 
enough to do the required matrix inversions. Despite 
these difficulties, recent developments from the field of 
numerical weather prediction 0, [f| 0, H, S [HJ E| sug- 
gest the possibility of achieving good accuracy (as in a 
Kalman filter), but in a way that is computationally fea- 
sible for large systems. 

The purpose of this Letter is to use numerical sim- 
ulations to demonstrate the potential usefulness of a 
new, weather-inspired data assimilation method, the lo- 
cal ensemble transform Kalman filter (LETKF), to a 
broad class of high dimensional spatiotemporally chaotic 
physical processes. For specificity we employ a par- 
ticular paradigmatic example: spatiotemporally chaotic 
Rayleigh-Benard convection. Flows such as spiral de- 
fect chaos [HI, [l3[ in the Rayleigh-Benard problem are, 
perhaps, the best studied experimental examples of spa- 
tiotemporal chaos; nevertheless, many general aspects of 
spatiotemporal chaos remain poorly understood. The 
LETKF is motivated by the observation that, in typ- 
ical examples of spatiotemporal chaos, spatial regions 



much smaller than the system size are accurately de- 
scribed by many fewer degrees of freedom than the full 
system. With this in mind, the LETKF employs many 
independent data assimilations in a large set of heavily 
overlapping regions. Because these regions are relatively 
small, the individual regional computations are not pro- 
hibitive. In addition, the regional data assimilation com- 
putations are independent of each other and can thus be 
done in parallel. Furthermore, by use of a simple exam- 
ple [3, [n| it was indicated that, if the size of the indi- 
vidual regions employed in LETKF is not too small (but 
still small compared to the total system size) , then state 
estimates with accuracies virtually the same as those for 
a classical Kalman filter technique (thus presumably of 
near optimal accuracy) can be achieved. For details of 
the LETKF algorithm we refer the interested reader to 
Refs. HESIll. 

In Rayleigh-Benard convection, a horizontal fluid layer 
of thickness d is confined between a heated lower plate 
and a cooled upper plate. The onset of fluid mo- 
tion occurs when buoyancy overcomes viscous dissipa- 
tion and thermal diffusion as the temperature differ- 
ence between the plates AT is raised above a critical 
value AT C . Rayleigh-Benard convection is typically mod- 
eled using the Boussinesq equations [13] , which are com- 
monly nondimcnsionalizcd with temperature scaled by 
AT, length scaled by d, and time scaled by the vertical 
diffusion time d /k, where k is the thermal diffusivity. 
This system of units is used throughout this Letter. We 
numerically solve the Boussinesq equations [1 81 ] applying 
realistic boundary conditions u = and T = (conduct- 
ing) on the walls of the region boundary (x 2 + y 2 ^ T 2 , 
\z\ ^ |). We denote the temperature deviation from the 
conducting static solution (which is linear in z) as T and 
the fluid velocity as u. 

The Boussinesq equations have two dimensionless pa- 
rameters, the reduced Rayleigh number e and the Prandtl 
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Here R = gacPAT/isn is the Rayleigh number, R c is the 
critical Rayleigh number at convective onset, g is gravi- 
tational acceleration, a is the thermal expansion coeffi- 
cient, and v is the kinematic viscosity. Fluid convection 
arises when e > 0. In addition, the radius T of the disk 
in units of the cell depth d, also referred to as the aspect 
ratio, specifies the geometry. We focused our studies on 
r = 20, e = l, Pr = l H. 

In experiments, flows are visualized using the shadow- 
graph method [20J, an indirect measurement of the fluid's 
spatial temperature variation. In typical experiments, 
due to its difficulty, measurement of the fluid velocity 
field is not performed. The so-called mean flow, defined 
here as u(x,y) = J uj_(x,y, z)dz (though see [2l| for 
a more complete description), has been shown, through 
the use of simulations, to play a significant role in the 
dynamics [21(. Here uj_ denotes the horizontal compo- 
nent of the fluid velocity u = u± + u z z. Our goal is to 
determine the full fluid state (Tlx, y, z), u(x, y, z)), from 
a time series of shadowgraph measurements at a finite 
number of horizontal (pixel) locations, and we view this 
as a test case investigation of the general usefulness of our 
technique for laboratory experiments on spatiotemporal 
chaos. 

We begin by considering a system state vector £ with 
N components, for which we have a dynamical model, 
= G(£j). Our G(-) is an integration of the Boussi- 
nesq equations 18] from a time tj to tj+i where tj = jAt 
and t\,t2, . . ■ are the times at which state estimates are 
constructed (also the times at which a measurement is 
taken). The model state £ consists of the variables T 
and u defined on the grid points of the cylindrical mesh 
used by the model. 

We map the temperature field to the shadowgraph light 
intensity I(x, y) with a map M using a relation derived 
from geometric optics [2(], HH 
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Here, Vj_ = d 2 / dx 2 + d 2 / dy 2 is the horizontal Laplacian, 
and the temperature field T(x, y) = J T(x, y, z)dz is ver- 
tically averaged. I (x,y) is the incident light intensity 
and a = 2zi\dn/dT\, where n is the index of refraction 
of the fluid, z\ is the optical path length from the mid- 
plane of the fluid layer to the image plane (in units of 
d) , and dn / dT is evaluated at the average fluid temper- 
ature. Note that for © to be a good approximation 
to a more correct physical optics treatment [23j we re- 
quire ||aV5_T|| 1, thus Af [■] is only weakly nonlinear in 
T(x,y,z). We assume that shadowgraph measurements 
are of the form M[T] mn + e mn where M[T] mn denotes 
the value of M[T(x, y, z)\ at the location of pixel (m, n), 



and the quantities e mn denote the errors in the measure- 
ment (noise), which we take to be zero mean iid Gaussian 
random variables with variance a 2 . 

Most data assimilation algorithms are iterative, cycling 
between a predict and update step once every time in- 
terval At. In the update step, current measurements are 
used to update (or correct) the prediction. The predict 
step then propagates the updated state, via the model, 
to the next measurement time (i.e., it is a short term 
forecast). If the method works as intended, this process 
will closely synchronize the experiment and the model 
by coupling them via the measurements. The LETKF 
operates on this basic principle [3, 0j| [l6| . 

In order to assess how well the LETKF is performing, 
we will compare it to a more naive approach that we call 
Direct Insertion (DI). At the time tj of the shadowgraph 
measurement Ij(x,y), the DI method updates the pre- 
dicted temperature field T^(x,y, z) by adding a correc- 
tion STj(x, y, z) that is the unique field that is quadratic 
in z, matches the boundary conditions at \z\ — i, and 
for which the updated field T™(x, y,z) = Tj(x,y,z) + 
STj(x, y, z) satisfies M[T"(x, y, z)] = Ij(x, y). This gives 
the update 

STjix^z) = (T?(x,y)-ff(x,y)) [(3/2) - 6z 2 ] , 
where T™(x,y) is found by solving a Poisson equation, 
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and (x c ,y c ) is the location of the closest pixel to (x,y) 
that is observed. Note that with DI the velocity is not 
updated, (u"(x,y,z) = u^(a;, y, z)); rather, it develops 
via coupling with the temperature during the simula- 
tion step. The predicted field (which has a proper z- 
dependence) is only slightly affected since, if the method 
is working properly and measurements are sufficiently fre- 
quent, the correction STj(x,y, z) is small. This method 
is the most successful data assimilation technique we 
have tested that does not use an update based on the 
Kalman filter. It is meant to represent what one might 
try without knowledge of the techniques described in 
Ref. 0, [H [3. 

Now we describe so-called perfect model tests in which 
a time series of states, generated from a Boussinesq simu- 
lation of one particular initial condition, serves as a proxy 
for the evolution of an experimental system whose state 
we wish to infer. Simulated shadowgraph measurements 
of this time series are generated every At — 1/4 using @ 
with the parameters a = 0.08, I Q (x,y) = 0.5, and adding 
noise of variance a 2 as previously described. Measure- 
ments are made sparse by removing shadowgraph pixels, 
leaving only those that lie on observation locations. We 
introduce the measurement density p = s/(irT 2 ), where 
s is the number of observation locations. When p is not 
small (p > 5) we randomly and uniformly distribute ob- 
servation locations over the disk; while at low density 
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(p < 5) the observation locations are placed on a Carte- 
sian grid covering the disk x 2 + y 2 ^ T 2 (giving more re- 
peatable results when using sparse measurements). The 
observation locations are fixed for the entire data assim- 
ilation process. 

We apply the LETKF and DI methods to our sim- 
ulated shadowgraphs to approximately reconstruct the 
original time series of true states. Here we document 
their performance as a function of measurement noise a 
and measurement density p. Simulated shadowgraphs are 
assimilated at times tj, j = 1 . . . J. During this process 
the DI and LETKF converge on an estimate of the system 
state (J chosen large enough to ensure convergence). At 
time tj assimilation is turned off and the final updated 
state estimate is used as an initial condition for a long 
term forecast. 

Performance is quantified via the RMS rela- 
tive error of the temperature and mean flow, 
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FIG. 1: The error of the forecast temperature Er(t) with 
a = 0.01 and p = 127. The dashed line is our chosen threshold 
defining the predictability time Eritp) = 0.15. 



1 

. 8 
. 6 



. 4 



0.2 




are the true temperature and mean flow from the proxy 
experiment, and (•) indicates a spatial average. The fields 
T(x, y, z, t) and u(a;, y, t) are the long term forecast fields £ n 
from either DI or LETKF. Note that, in contrast with 
the case of a real physical experiment, the perfect model 
set-up used here has the advantageous property that the 
exact true fields u* and T* are known and hence available 
for evaluating the actual error Et and E a . 

Wc define the ideal scenario as measuring a shadow- 
graph with p = 127 (corresponding to a 451 x 451 shad- 
owgraph image) and a = 0.01 (this situation can be 
achieved in an experiment). Under these conditions the 
DI and LETKF converge on a state estimate within a 
few vertical diffusion times. Both DI and the LETKF 
are effective for estimation of the (unobserved) mean flow 
u(x, y); however, the LETKF achieves an initial error E^ 
that is less than half that of DI. The forecast error for a 
typical state estimate is shown in Figs. [T] and O The gen- 
eral character of the forecasts is near-perfect agreement 
with the true state, followed by rapid divergence due to 
local error growth at the location of a defect. It is clear 
that the LETKF forecast is far superior to DI's, as mea- 
sured by the predictability time t p , defined as the time 
when the forecast error Et first crosses the (somewhat 
arbitrary) value of 0.15. 

Under non-ideal conditions the LETKF proves much 
more robust than DI. Results for sparse measurements, 
shown in Fig. [3J demonstrate the large range of p for 
which the LETKF converges. One can observe the exis- 
tence of a critical density of observations (p w 1.3) below 
which it fails to converge. DI on the other hand exhibits 
a rapidly deteriorating forecast when even a few obser- 
vation locations are removed. 

We investigated performance as measurement noise G(£j,Pj) , 

was increased. The meaningful signal to noise ratio is Tj+i — Pj+i Pj ~ 

<j S g/o~, where a sg = (I{x, y) — (I{x, y))) is the RMS inten- 
sity of a typical shadowgraph (a sg — 0.123 when a — 0.08 By now regarding the system model as G, estimates of 
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FIG. 2: The error of the forecast mean flow Ea(t) (an unob- 
served variable) with a = 0.01 and p — 127. 



and I (x,y) = 0.5). DI relies on the Poisson solve ([3]) 
which is fundamentally insensitive to noise (it smoothes 
the right hand side). However, this insensitivity com- 
petes with the sensitivity of the chaotic system dynamics 
when producing forecasts. The net result, in Fig. 2] indi- 
cates that DI forecasts are only reliable for a few vertical 
diffusion times when a > 0.4cr sg , whereas the LETKF 
yields accurate forecasts up to and exceeding a = a sg . 

Finally, we note that it is common for model param- 
eters to be unknown, and that one of the advantages 
of the Kalman filter methodology is that it can be uti- 
lized to infer unknown system parameters from measured 
time series [lfSj]. In particular, let p denote the vector of 
unknown system parameters; with the model now deter- 
mined by — G(£j,p). One can then introduce an 
extended state space vector having the form 7 = [£ p] T , 
where p is treated as a state variable with no time de- 
pendence. The extended model evolves as 
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FIG. 3: Comparison of DI and the LETKF; the predictability 
time t p is shown as the density of observations p is reduced, 
demonstrating the superiority of the LETKF. 




0.5 1 1.5 2 



FIG. 4: The predictability time t p is shown as measurement 
noise is increased. Noise levels shown are a = [0.01, 0.05, 0.08] 
for DI and a = [0.01, 0.05, 0.2] for the LETKF. 



7 (and therefore the parameters p) result from an im- 
plementation in the same way as for but in the space 
of 7 vectors. Using this method we can achieve esti- 
mates of e to within 0.02% of the true value (e = 1) 
with the LETKF. Remarkably, even when measurements 
are extremely sparse (e.g., p — 3.6, near the critical 
measurement density) the e estimates are within 0.2%. 
This demonstrates the utility of the LETKF for inferring 
model parameters from incomplete data. 



In conclusion, our results support the potential effec- 
tiveness of the LETKF at estimating the fluid state in 
laboratory Rayleigh-Bcnard convection experiments. We 
believe that the method we have presented is applicable 
to a large class of spatiotemporally chaotic systems, and 
offers the possibility of bridging gaps between experiment 
and theory in the study of spatiotemporal chaos [24[ ■ 
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